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A theoretical model of a single molecule coupled to many vibronic modes is presented. At low 
energies, transport is dominated by electron- vibron processes where transfer of an electron through 
the dot is accompanied by the excitation or emission of quanta (vibrons). Because the frequency 
of the nth mode is taken as an nth multiple of the frequency of the fundamental mode, several 
energetically degenerate or quasidegenerate vibronic configurations can contribute to transport. We 
investigate the consequences of strong electron- vibron coupling in a fully symmetric setup. Several 
striking features are predicted. In particular, a gate- asymmetry and pronounced negative differential 
conductance features are observed. We attribute these features to the presence of slow channels 
originating from the interplay of Franck-Condon suppression of transport channels with spin and/or 
orbital degeneracies. 

PACS numbers: 73.23.-b,85.85.+j,73.63.Kv 



I. INTRODUCTION 



The coupling of electronic and mechanical degrees of 
freedom is at the core of the physics of nanoelectrome- 
chanical systems (NEMS) 1 . Pronounced vibrational ef- 
fects in electronic transport have been observed in several 
recent experiments on molecules 2 - — and single wall car- 
bon nanotube quantum dots 7 - - — . Stimulated by the ex- 
perimental works, several groups— ~— have attempted to 
theoretically explain some features of the measured sta- 
bility diagrams (i.e., of the differential conductance in a 
bias voltage-gate voltage colour map) which appear to be 
ubiquitous. Specifically, one observes the following: (i) 
Equidistant lines run parallel to the edges of the Coulomb 
diamonds, (ii) The probed diamonds show often nega- 
tive differential conductance (NDC) features appearing 
between excited vibronic states, (iii) This corresponds 
to a sequence of peaks in the current as a function of 
the bias voltage. The height of the peaks measured in 
Ref. was found to be in quantitative agreement with 
predictions of a simple Franck-Condon model for a sin- 
gle electronic level coupled to a harmonic mode (the so 
called Anderson-Holstein model^^i 7 .. The NDC fea- 
tures have been explained in Ref. [3 in terms of local 
vibronic excitations. In Ref. UH NDC features are as- 
sociated with the Anderson-Holstein model at a single 
vibronic resonance, while in Ref. LL9| NDCs appear due 
to the electron-phonon-coupling-induced selective unidi- 
rectional cascades of single-electron transitions. In a re- 
cent work—, carbon nanotube-specific NDC features have 
been attributed to a spatially dependent Franck-Condon 
factor and to the presence of quasidegenerate electronic 
levels, while in Ref. [lO an interference effect between the 
orbitally de gen erate electronic states has been proposed. 
In Refs. [I2I.T19I. and [lO, however, to see NDC some sort 
of asymmetry is required either in the coupling to the 
source and drain contacts or in the coupling to the two 
orbitally quasidegenerate states of the system or in both 
couplings at the same time. 

In this paper we extend these ideas and propose a 



generic model in which two degenerate or quasidegener- 
ate molecular levels are coupled to many vibronic modes. 
We consider the case of a symmetric setup, i.e., invariant 
under exchange of source with drain if simultaneously 
also the sign of the bias voltage is reversed. If the fre- 
quency uj n of the nth mode is an nth multiple, uj n = nou, 
of the frequency uo = uo\ of the fundamental mode, then 
several energetic degenerate vibronic configurations arise 
(involving two or more modes) which may contribute to 
transport at finite bias. 

As we are going to show in our paper, the additional pres- 
ence of spin and/or orbital degeneracies opens the possi- 
bility of getting slow channels contributing to transport. 
As a consequence, NDC phenomena can occur despite a 
fully symmetric quantum-dot setup. A peculiarity of the 
observed features is, in particular, an asymmetry with 
respect to the gate voltage in the stability diagrams. Fi- 
nally, we retain source-drain symmetry but allow for an 
asymmetry in the coupling to the two orbitally degener- 
ate states. We show that this asymmetry is sufficient to 
explain the experimental observations (i)-(iii). 

In some of the experiments 7,12 the slope of the NDC 
lines is the same for both positive and negative biases. 
This characteristic has been associated with the left and 
right asymmetry in the coupling to the leadsi 2 -. We con- 
firm that including the higher harmonics does not change 
the effect and give an analytical interpretation of the nu- 
merical results for low biases. 

The paper is organized as follows: In Sec. II the model 
Hamiltonian of a single molecule coupled to several vi- 
bronic modes is introduced. A polaron transformation 
is employed to obtain the spectrum of the system in 
the presence of electron-vibron interactions. In turn, as 
known from the theory of Franck-Condon blockade in 
the simplest Anderson-Holstein model 1 the polaron 
transformation also crucially affects the tunneling Hamil- 
tonian describing the coupling to the source and drain 
leads. 

In Sec. Ill the consequences on transport are analyzed, 
(i) The tunneling transition amplitudes involve product 
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of Franck- Condon factors with coupling constants de- 
pending on the mode number, (ii) At low bias, such that 
only the lowest vibronic mode is excited, a description of 
the dynamics only in terms of rate equations involving 
occupation probabilities of the many-body states of the 
quantum dots is appropriate, (iii) At higher bias, when 
several vibron modes are excited, a generalized master 
equation (GME) coupling diagonal (populations) and off- 
diagonal (coherences) elements of the quantum dot re- 
duced density matrix should be used (see e.g., Refs. I20I - 

In Sec. IV our main results for the current- volt age 
characteristics in a fully symmetric setup are shown and 
analyzed. In particular, we give an explanation of the 
NDC features observed at low bias in terms of a differ- 
ent spin and orbital degeneracy of states with different 
electron numbers contributing to transport. In Sec. V 
asymmetric setups are discussed as well. Conclusions are 
drawn in Sec. VI. 



II. MODEL HAMILTONIAN 

In the following we consider a generalized Anderson- 
Holstein model where the Hamiltonian of the central sys- 
tem has the form 



H, 



sys 



mol 



Hy + H e 



(1) 



where Hmol describes two quasidegenerate levels. The 
Hamiltonian is modeled as 

Hmol = + ^ - X ) ' ( 2 ) 

la 

where / = 1,2 is the orbital and a =t, 4- is the spin 
degree of freedom. The operator Ni a = dl a dia counts 
the number of electrons with spin a in the orbital /. 
N = ^2i a Nia is the total number operator. The or- 



bital energy is e\ = So 



with A an orbital- 



1 + (-1) A 

mismatch. The Coulomb blockade is taken into account 
via the charging energy U and we assume U > £q. 
The vibron Hamiltonian is expressed as 

V 



a'a ri 



(3) 



where a n (ajj annihilates (creates) a vibron in the nth 
mode of energy e n = hw n . We assume that the energy of 
the nth mode is given by 



nhu, 



(4) 



being an nth multiple of the energy e\ = Jtlj of the 
fundamental mode as it is, for example, for longitudi- 
nal stretching modes in quantum wires and carbon nan- 
otubes. Finally, the electron- vibron interaction Hamilto- 
nian is given by 



H e - V = ^2^2gnNi a (d ] n 



(5) 



n>i No- 



where g n is the coupling constant for the nth vibronic 
mode. 



A. Polaron transformation 

In order to solve the Hamiltonian of the system, we 
decouple the electron- vibron interaction part by applying 
a unitary polaron transformation^. Specifically, we set 

H sys = e s H sys e~ s , where 



S =^2^2\ n N la (fit - On) 



(6) 



n>l la 



and A n = g n /hw n is the dimensionless coupling constant 
associated with mode n. Notice that A = is the 
coupling constant of the fundamental mode. We assume 
A = 0.68, 0.83 and 1.18 in the analysis of the spec- 
trum, which is in the range of values observed e.g., in 
experiments on carbon nanotubes^^ii^. Under the po- 
laron transformation the operator d a \ is transformed as 



di a = e s di a e~ 



di a X, 



(7) 



where X = exp 



En>l ^ n fan ®n 



. In a similar way 



the shifted vibronic operator is 

d n = d n — \ n Nia- 



(8) 



lo 



Ha 



The transformed form of the system Hamiltonian is thus 

la n>l ^ 

^n(n-i), (9) 



where si = si — ^2 n JjN- is the renormalized orbital en- 

~ n i i 2 

ergy and U = U — 2 j^j- is the Coulomb repulsion 

modified by the vibron mediated interaction. 

The eigenstates of the system are 



\N,m v ) 1 :=e s \N,m v ). 



(10) 



where N = (Ni^, N±i, N2^, N21) and Ni a the number of 
electrons in the branch (la). Notice that N = Ni a 
defines the total number of electrons on the dot. For later 
purposes we indicate the ground state and first excited 
state with electrons as [see Fig. EJb)] 



10,0} := |0,0>!, 

|0,1) := |0,m„ = (l,0,0,...))i, 



(11) 



The first excited state with 0-electron contains one vi- 
bronic excitation in the first mode, i.e., m v = (1, 0, 0, ...). 
In a similar way we define the ground states and first 
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excited states with N = 1 electron. For zero orbital mis- 
match one has fourfold degeneracy [see again Fig. E^b)], 
i.e., 



|l fc ,0):=|ffc,3)i, k = 1,2,3,4, 
|l fc ,l) := \U,m v = (1,0,0,...))!, 



(12) 



where l k e {(1, 0, 0, 0) , (0, 1, 0, 0),(0, 0, 1, 0) , (0, 0, 0, 1)}. 
We notice that both for N = and N = 1 the second ex- 
cited states are vibronically degenerate for the dispersion 
relation Eq. ((4]), then the configurations rh v = (2, 0, 0, ...) 
and m v = (0, 1, 0, ...) have the same energy. For finite or- 
bital mismatch, the orbital degeneracy is broken. 
The case < eoA < Huj is illustrated in Fig. EJb). The 
corresponding states are: 



llfcfl.0) 

|lfce,0> 

llfcfl.l) 
|lfce, 1) 



|l feg ,0)i, k = l,2, 
|l fee ,0)i, k = l,2, 
\lk g ,m v = (1,0,0, 
\Ue,m v = (1,0,0, 



(13) 



and l kg e {(1,0,0,0), (0,1,0,0)}, 
f fce €{(0,0,1,0), (0,0, 0,1)}. 



degrees of freedom is taken. We make the following 
standard approximations to solve the above equation: (i) 
The leads are considered as reservoirs of noninteracting 
electrons in thermal equilibrium. (ii) We factorize 
the total density matrix as p 1 (t) = pi ys (i) <S> pleads, 
where we assume weak coupling to the leads, and treat 
Ht perturbatively up to second order. (iii) Being 
interested in long time properties, we make a Markov 
approximation, where the time evolution of f>l ed (t) is 
only local in time. Notice that in the regime of interest, 
t — > oo, the Markovian approximation becomes exact, 
(iv) Since the eigenstates \N,m v } 1 of H sys are known, 
it is convenient to calculate the time evolution of pl ed 
in this basis retaining coherences between degenerate 
states with the same number of particles. Hence p T red 
can be divided into block matrices Pn^ N (i), where 
E and N are the energy and number of particles of 

the degenerate eigenstates \n), \m) G ||7V, m v ) 1 |. We 
obtain an equation of the Bloch-Redfield form 



Vm^W — ^ Rnmkk' Pkk' ** ft) 



£ ( 17 ) 



M=7V±1 E' kk' 



III. SEQUENTIAL TRANSPORT 

In this section we discuss the transport across the sys- 
tem under the assumption of weak coupling to the leads. 
The Hamiltonian of the full system is described as 



H — H s 



Hfy + Hj 



Hp 



(14) 



where a = s,d denote the source and the drain contact, 
respectively. The tunneling Hamiltonian Ht is given by 



olkct I 



H.c.) 



(15) 



where c aKa is the electron operator in the leads. Finally, 
H e xt describes the influence of the externally applied 
gate voltage V g . The gate is capacitively coupled to the 
molecule and hence contributes via a term eV g N. In the 
case of high degeneracy of the spectrum, the appropriate 
technique to treat the dynamics of the system in the weak 
coupling regime is the Liouville equation method for the 
time evolution of the density matrix of the total system 
consisting of the leads and the generic quantum dot. To 
describe the electronic transport through the molecule, 
we solve the Liouville equation 



ih 



dt 



Tvi ea d s 



'Ufa), fit) 



(16) 



for the reduced density matrix p re d(t) — ^ T ieads {pft)l 
in the interaction picture, where the trace over the leads 



where the indices n, m, fc, k' refer to the eigenstates of 
H S ys and fc, k' runs over all degenerate states with fixed 
particle number. Notice that if TV = then M = 1, while 
if N = 4 then M = 3. The Redfleld tensors are given 
by^ 



B E N _ V- y> (g T (+)E N E' M g T (-)E N E' M \ 
■ — ikk' ~ y u mk"- a,njjk ' u nk L a ^k' jjm )i 



a=s,d M,E',j 



(18) 



and Rllp = E Q , P=± rS?if^ where the quantities 

T^njjk™ are transition rates from a state with N to a 
state with M particles. More explicitly: 



(p)e n e' n+1 _ y- r (p) , ) (d, \ EnE ' n+1 ($ Y 

la 



E' N+1 E N 



(19) 



with e a = — eV a — (En — ^yv+i) an d Va the electrochem- 
ical potential of the lead a. Likewise 



(p)E N E' N _ 1 
a,k' rank 



EnE' n _ 1 / „ \E , n _ 1 En 



Erg-(4)(4)r-'W 

la 



(20) 



with e' a = — eV a — (E' N _ 1 — E^. Moreover, we intro- 
duced 



rJL (E) = lal f± (E) + l ^ lal P J de 



f±{e) 
e-E' 



(21) 



where /+ (e) = f (e) is the Fermi function while /_ (e) = 
1 — / (e) and j a i = jfD a \t a i\ 2 are the bare transfer 
rates with the constant densities of states of the leads D a . 
Knowing the stationary density matrix pl t , the (particle) 
current through lead a is determined by (a = s/d) 

j _ V" (t WEnE ' n + 1 - rW^^-^ J, En 

l a - ^a/ie / J \ L ot.njjk 1 a,njjk J Hkn,at' 

N,E,E' nkj 

(22) 

If the relation given by Eq. (j4j) holds, then spin and 
orbital degeneracies intrinsic in the electronic structure 
are supplemented by degeneracies related to the vibronic 
structure. Several vibronic modes with frequencies uo n = 
nuj multiples of the fundamental frequency uo give rise, in 
fact, naturally to several degenerate vibronic configura- 
tions. This is the situation we shall focus on in the rest 
of the paper. 

A degenerate spectrum is a necessary condition for the 
appearance of interference effects in the transport char- 
acteristics both in the linear and non linear regime^ — 
and these effects can be captured only by considering not 
only populations (diagonal elements) but also coherences 
(off-diagonal elements) of the reduced density matrix. 

For the system at hand we calculated the current both 
with and without coherences between degenerate states 
up to five vibronic modes, obtaining though only quan- 
titative but not qualitative differences. While spin and 
orbital degeneracies can be a priori excluded from the 
transport through a single molecule with nonpolarized 
leads^i, the role played by the vibronic coherences re- 
quires a more careful analysis. 

We have confirmed that it is not possible to construct 
a linear combination of degenerate states {\s}} with fi- 
nite transition amplitude to a state |r) at one lead but 
decoupled from \r) at the other lead, where \r) and \s) 
represent the states given by Eq. (fTD|). This observa- 
tion, complemented by the general method presented in 
Ref. [33| proves the absence of interference blocking states 
in our system. Thus, interference, even if present, does 
not have dramatic consequences on the transport char- 
acteristics of the system. 

All the current maps presented in the next section are 
hence, apart from Fig. [61(b), calculated neglecting coher- 
ences. As shown explicitly in Fig. [6fb), this approxima- 
tion does not affect qualitatively the results (at least in 
the low bias regime). Moreover, the negative differen- 
tial conductance and the associated dynamical symme- 
try breaking that we present in the next section are not 
related to the interference and can thus be obtained by 
considering the dynamics of the populations alone. 

IV. SYMMETRIC SETUP 

In this section, we illustrate our predictions for the 
transport characteristics and focus on the <H> 1 transi- 
tions. In the calculation we also assume for the coupling 



constant of the nth mode g n = \/ngi (as expected for 
stretching modes in carbon nanotubes^i) . The system 
is symmetrically coupled to source and drain contacts 
(7sZ = Idi) and the lowest five vibron modes are included. 
Results for the differential conductance for different val- 
ues of the dimensionless coupling constant A are illus- 
trated in Fig. Q] and Fig. [7J obtained for zero and finite 
orbital mismatch ea = £oA, respectively. 

A. I — V characteristics at low bias and zero band 
mismatch 

When the orbital mismatch is zero, i.e., sa = 0, then 
the two orbital energies e\ are the same. The minimum 
energy to produce a charge excitation is eq. We take the 
value of this energy as So = lAmeV (comparable to the 
level spacing energy of a suspended single wall carbon 
nanotube of 1.2/ira length). Furthermore, we assume the 
energy of the lowest vibronic mode to be E\ = 0.04meV. 
Thus the charge excitation energy is much larger than 
the energy of the lowest vibron mode. Indeed all the 
equidistant lines running parallel to the diamond edges 
observed in Fig. [T] are due to vibron excited states. What 
striking is the occurrence of negative differential conduc- 
tance (NDC) features at moderate coupling (A = 0.68 
and A = 0.83) which, however, disappear when the cou- 
pling is increased (A = 1.18). Moreover, the NDC lines 
are only running parallel to one of the diamond edges, 
which indicates an asymmetry with respect to the gate 
voltage V g . As we are going to explain, at low bias, 
these features are a consequence of Franck- Condon as- 
sisted tunneling combined with the spin and/or orbital 
degeneracy in the system. 




eV g [meV] eV g [meV] ' eV g [meV] 

Figure 1: (Color online), (a)-(c) Plots of the numerical dif- 
ferential conductance dI/dV( arbitrary units) of the system 
for coupling constants A = 0.68, 0.83 and 1.18, respectively. 
The charge excitation energy is so = 1.4meV and the energy 
of the lowest vibron mode is si = 0.04meV. Additional pa- 
rameters are a thermal energy of ksT = 0.8/xeV, orbital mis- 
match sa = and ^ s i = 7dz = 0.02/xeV for / = 1, 2. The black 
lines running parallel to the Coulomb diamond edges cor- 
respond to negative differential conductance (NDC). Notice 
that here and in the following figures dI/dV( arbitrary units) 
is normalized to the maximum of dl/dV (arbitrary units) in 
the considered parameter range. The gate voltage is set to 
zero by convention at the degeneracy point. 

Specifically, let us focus on the low bias region [see 
Fig. [51(a)], where only ground-state <H> ground-state tran- 
sitions (region A), and ground-state o first excited-state 
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transitions (regions B, C) are relevant. The and 1- 
particle states involved are illustrated in Fig. Efb), to- 
gether with their degeneracy due to spin and orbital de- 
grees of freedom, and have energies below the dashed 
line in Fig. [2fb). The states above the dashed line re- 
quire an energy of at least 2hu and have thus also a 
vibron degeneracy. In the considered energy range no 



— (2x) — (8x) 




_ (ix) (4x) 

:(1x) -rrr(4x) 



N=0 



N=1 



and the four 1-particle ground states |lfe, 0) [see Eq. ([12]) ] 
with energy E® contribute to transport. Moreover, inside 
region A it holds 



f(eV s -(E°-E° )) = l, 
l-f(eV d -(E° 1 -E° ))=l, 



(26) 



such that, if 7 a i = j a2 = 7a, 7s 
(|n> = |0,0), |fc>e{|l fc ,0») 



7d, it also follows 



a,nk / j a,kn 1 



00- 



(27) 



This situation is illustrated in the table of Fig. [3j where a 
dashed red (black) arrow indicates a transition involving 
the source (drain). Condition ([27]) with ([24]) then implies 
that 



Figure 2: (Color online), (a) The low-bias transition regions 
of the stability diagram are labelled as A, B, C, D. (b) Energy 
level scheme for the relevant transitions in the stability dia- 
gram involving regions A-D. Above the dashed line region D is 
activated with two vibron modes being in the transport win- 
dow. The energy of the lowest vibron mode is hw — 0.04meV. 
The number of degenerate states is indicated in bracket. 

degenerate vibron configurations are involved. Moreover, 
coherences between degenerate electronic configurations 
are not present such that a rate equation description only 
in terms of populations is appropriate. At low bias and 
in the stationary limit Eq. ([T7]) yields the equation for 
the populations: 







p En 

rnn 



-R. 
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E EE*' 



M=7V±1 E' 



ikk rkk ' 

(23) 



E E 

pnn = p k k Vfe, and n 



(28) 



and hence Pfi '.= pnn = \] Pf = Y,k Pkk = f> yielding 
with Eq. ([25]) for the current in region A, I a = fToo- 
Along similar lines we can calculate the current in re- 
gions B and C. Let us start with region B where (see 
table of Fig. [3]) the gate voltage V g is such that the 
1-particle ground states |lfe,0) have energy E± smaller 
than the one, of the 0-particle ground state |0,0). 
Moreover, in this region also the first excited state |lfc, 1) 
with energy E\ enters the transport window. We also 
assume that the rate Tn between the states |0, 1) and 
| lit, 1) is negligible with respect to Too and Toi. Correc- 
tions due to a finite Tn will be discussed later. Inside 
region B it holds, besides Eq. ([26]) . and hence Eq. ([28]) . 
/ (eV s - (El - E° )) = 1 and 1 - / (eV d - {E\ - E° )) = 
1. Hence it follows that (\n) = |0,0), \k) e {|l fc ,l)}) 



or, equivalent ly, 



E EEE(#^- r 



Pkk 



M=N±1 E' 



where T En J' m 



(24) 



ZUeL a,nkkn anCL 1 a, fen 



2ReT{ alnnk EN i see E( f s - ©-(ED- Notice in partic- 
laif+ (eV a — (E f N+1 — En)) C n u 



ular that rf^ 1 

and T E ^ En = lal f. (eV a - (E' N+1 - E N )) C nk , i.e., 
they only differ in the Fermi factors. The transi- 
tion coefficients are the same and given by C n k = 

^2 la ( dia ) • Moreover, we find for the current 

V /fen 

through lead a 

J «=* E EE(C? +1 - r Sr" w )^- (25) 



N.E.E' 



1 01 = / A cx.nk -Z^ L a,kn = 1 10- 



Eq. ([24]) implies thus that in region B it holds 



^0 

Pnn 



E» 
Pkk 



Pkk ~ g> 



(29) 



(30) 



and hence Pq = 



^0 

Pnn — 



E° 

Sfe Pkk 



EE 1 
k Pkk 



|. The total current in region B follows from 



Eq. ([25]) and reads Is = 9 (Too + Toi). The condition to 
observe NDC is that Ib < Ia, which implies 



Toi < ^roo- 

5 



(31) 



Along similar lines (see table in Fig. [3]), one finds for the 
transition from region A to C that Iq < I a if 



Toi < ^Too- 
5 



(32) 



Let us now focus on region A. In this case only the 0- 
particle ground state |0,0) [see Eq. ([TT]) ] with energy E® 



Let us look in more detail at Eq. ([31]) and Eq. ([32]) . 

The rates Too and Tqi describe transitions between states 
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which only differ in their vibronic part. From Eqs. ([27]) , 
and (|Aip it follows that 



^ = F 2 (A,0,1) 
J- oo 



A . 



(33) 



Hence, to observe NDC for the transition from region A 
to B one needs that A 2 < |. On the other hand for NDC 
in the transition from A to C we must require A 2 < |. 
Indeed, as shown in Fig. [TJ NDC for the transition A <H> 
B is observed for A = 0.68 and A = 0.83, but it vanishes 
for A = 1.18. On the other hand NDC is never observed 
for the transition region A «->> C. Let us now turn to 
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Figure 3: (Color online). In the table for each of the differ- 
ent regions A, B, C in the stability diagram, relevant tran- 
sitions, population of states and current are given. Pq , P e 
represent the population of the 0-particle ground and first 
excited states, respectively, and Pf, Pf the population of 
the 1-particle ground and first excited states. / is the cor- 
responding current in each region. In the transition scheme, 
the black arrows represent the drain and the dashed red ar- 
rows the source transitions. Too denotes the transition rate 
from 0-particle ground state to the 1-particle ground state 
while Toi the transition rate from the 0-particle ground state 
to a 1-particle first excited state. 



1^1 Er, 



region B and to a finite T n ee E« T Z°nk = T a K]L 
with \n) = |0, 1), k G {|lfe,l)}. Because now |0, 1) can 
get populated, also transitions from |0, 1) to 1 1/^,0) are 
activated (see Fig. 2]). 

— Ei > eV s ,eV d it holds 



Because of Eq 



10 



^01 = ^2*5$ = °- 



Hence, the stationary solution with equal probabilities is 
spoiled, at finite Tn, due to the inequality of Toi and 
fio- In fact, in the case Tn = 0, the same inequality 
only implies that p E o = 0. We also notice that Tiq = 
Tiq. Moreover, the rates Tn and Too only differ in their 
vibronic configuration: it holds [cf. Eq. (|A2[) ] 

2 

^ 2 



£11 
Too 



£(-a 2 )' 



Likewise 



En 
r i 



(1-A 2 V 
A 2 



Hence, if |A| « 1 it is indeed Tn <C rocToi and an 
expansion to lowest order in the ratios ^P- can be 
performed. In this case the conditions for acquire 




Figure 4: (Color online), (a) Energy level scheme for transi- 
tions in region B and (b) in region C. Importantly, because 
the bias voltage is too low, the transition |lfc,0) — »• |0, 1) in 
region B and the transition |0,0) — > |lfc,l) in region C are 
not allowed. 

a more complicated form. The condition to get NDC in 
the source threshold lines [from A to B in Fig. [2a)] is 



Foi< 5 r °°"40 rn ' 



(34) 



while the condition for NDC in the drain threshold lines 
[from A to C in Fig. is 



Toi < ^Too 
5 



7_. 

To 



(35) 



It means that the presence of chain transition processes 
redistributes the population among the many-body states 
in a way that privileges the low energy states (see Fig.[5j); 
this in turn weakens NDC since it privileges the conduct- 
ing channels that carry more current. Eventually, let us 
consider explicitly the effects of the higher harmonics and 
of the coherences between states with different vibronic 
configuration on the transport characteristics of the sys- 
tem. In Fig. [6ja) we present the stability diagram for 
a coupling constant A = 0.68 in which we artificially 
neglect the higher harmonics. By a direct comparison 
with figure Fig. QJa), it is clear that this approximation 
only marginally affects the NDC and positive differential 
conductance (PDC) pattern, thus confirming the dom- 
inant role played by the spin and pseudospin (orbital) 
degeneracies in the gate asymmetry. The effect of the 
coherences, shown in FigE^b), is more complex. Noth- 
ing changes for the lowest transition lines where no de- 
generacy is involved. For higher biases, though, some 
drain transition lines change their character from PDC 
to NDC. Thus, the gate asymmetry introduced by the 
spin and orbital degeneracy and the corresponding NDC 
(PDC) character of the source (drain) transition lines is 
exact in the low bias limit but should be taken only as a 
trend when several excited vibronic states participate in 
the transport. 



B. I — V characteristics at low bias and finite band 
mismatch 

In this section, we discuss our results on vibration- 
assisted transport with the same parameters as in 
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Figure 5: (Color online). Populations of the low energy states 
corresponding to the stationary density matrix calculated for 
different electron vibron coupling A and different gate-bias 
ranges. The first column corresponds to the case A = 1, thus 
Tn = 0, while in the second column A = 0.83. The letters A, 
B, C labeling the rows refer to the stability diagram regions 
defined in Fig. [2] The states are ordered in energy. The rest 
of the parameters are the same as used for Fig. [1] 
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Figure 6: (Color online). Plots of the differential conduc- 
tance for a coupling constant A = 0.68 with two different 
approximations: (a) neglecting the higher harmonics of the 
system vibrations, (b) keeping coherences between the degen- 
erate states with different vibronic configurations. The rest 
of the parameters are the same as used for Fig. [Ha). 
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Figure 7: (Color online), (a)-(c) Stability diagrams for cou- 
pling constants A = 0.68, 0.83 and 1.18, respectively. Addi- 
tional parameters are a thermal energy of k&T = 0.8/xeV, or- 
bital mismatch sa — 0. 016eo = 0.56^i and 7 S = 7^ = 0.02/xeV 
while for (d)-(f) sa = 0. 006^0 • The rest of the parameters 
are the same as used for Fig. [T] 



Sec. II V Al but with a finite orbital mismatch, i.e, ea ^ 0. 
In this case the orbital degeneracy is broken. The cor- 
responding stability diagrams are shown in Fig. [7J The 
analysis for the NDC conditions at low bias remains al- 
most the same as before. Slight differences occur because 
in this case the orbital degeneracy is lost and the popu- 
lations are redistributed over the many-body states in a 
different way. In Fig.[8{a), the different transition regions 
of the stability diagram have been labeled while the en- 
ergy level scheme has been shown in Fig. |5Jb), where the 
degeneracy of each state is given in brackets. We again 
truncate the process at the dashed line to analyze the 
lowest-energy excitations. As before I = 1 and I = 2 are 
the orbital degrees of freedom. The transition scheme 
for regions D and F is shown explicitly in Fig. [9] The 
current-voltage characteristics Fig. Efa)-(c) and [Tfd)-(f) 
are qualitatively the same as far as the mismatch is in the 
moderate regime k&T <C £a < ^ 5 the only difference be- 
ing the position of the resonance lines, that depends on 



the specific position of the energy levels. In other terms, 
despite the size and the position of the regions of the sta- 
bility diagram (Fig. [7j) depend on the mismatch £a, the 
value of the current in each region is independent of it. 
Thus a unified treatment of the two cases presented in 
Fig. [3 is allowed, despite their apparent qualitative dif- 
ferences. In particular, we can observe that the current 
in region B is larger than the current in region A since 
a new transport channel is opening with the same geo- 
metrical coupling (Too) when passing from A to B. This 
implies that the first source threshold transition line is al- 
ways a positive differential conductance (PDC) line. The 
current in C is equal to the one in A since, due to en- 
ergy conservation, no new transport channel is opening 
when passing from A to C. The corresponding resonance 
line is thus invisible in the stability diagram (see Fig. [7]). 
The transition that defines the threshold line separating 
A from C (| l^ e , 0) o |0, 1) at the drain) involves, in fact, 
states which are not populated in that bias and gate volt- 
age range. Finally, the comparison between currents in 
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Figure 8: (Color online), a) The low-bias transition regions 
are labelled as A, B, C, D, E, F, G. b) Energy level scheme for 
the transitions relevant in the low bias regions of the stability 
diagram. As in Fig. [2] is fiw = 0.04meV. The degeneracy of 
each state has been shown in brackets. 



the adjacent B and D regions and between the currents 
in the C and F regions results in conditions for the ap- 
pearance of NDC lines which are very similar to the one 
in absence of mismatch. In particular, the condition for 
NDC at the transition between regions B and D is iden- 
tical to the one for the transition between regions A and 
the B corresponding to zero mismatch given in Eq. 
The NDC condition for the transition between region C 
and F reads instead 



r i < 



V57-1 
28 



oo 



rn, (36) 



to be compared with the one for the transition between 
the regions A and C and zero mismatch given in Eq. (l35|) . 
A similar analysis can be repeated for higher-energy tran- 
sitions which participate in the transport for higher bi- 
ases. It is already clear though from the low energy 
transitions that a moderate breaking of the orbital de- 
generacy introduced by the finite mismatch e& does not 
change qualitatively the transport characteristics of the 
system. In particular, it preserves the presence or ab- 
sence of asymmetric NDC lines as a function of the elec- 
tron vibron coupling A (compare Figs. [Hand EJ). 




Figure 9: (Color online), (a) Transition scheme for region 
D at sa — 0.016^0 • (b) Transition scheme for region F at 

£A = 0. 016^0- 



asymmetry is introduced in the tunneling coupling of the 
molecule dot to the source and drain (75,1/2 = 7^,1/2)- 
The theory can produce, though, alternating PDC and 
NDC traces as discussed in Ref. [H if we assume that 
coupling of the I = 1 and I = 2 orbitals to be different. 
In Fig. [10] we have plotted the differential conductance 
(dl/dV) for an asymmetry parameter a = = 1/45 
where "1" , "2" represent the orbital degrees of freedom, 
respectively and a means source or drain. For conve- 
nience, in the numerical calculations, we use parameters 
as in Figs. [1] and [71 As seen by comparing Fig. [7] with 
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Figure 10: (Color online), (a)-(c) Stability diagrams for a 
molecule for the case of a coupling to the leads which depends 
on the orbital degree of freedom. All the parameters are the 
same as used in Figs. [T]and[3 The asymmetry is a = 1/45 
with orbital mismatch sa — 0.016so- 

Fig- OH at A = 1.18 NDC can now occur. Moreover, an 
alternation of PDC with NDC lines, as seen in the ex- 
periments^^ 2 - occurs. Repeating the same analysis as in 
Sec. UVBl we indeed find that: (i) the transition from re- 
gion A to B gives a NDC line for a > 3/2 independent of 
the value of A. (ii) the condition governing the transition 
from region B to D is now modified to be (at Fu =0) 



F (2) 
1 01 
P (2) 
1 00 



X 2 < 



2a 



5a 



(37) 



which, for a < 1, increases the range of A giving NDC 
also for A = 1.18 [Fig. HDTc)]. r$ and are defined 
analogously Too and Toi by considering the I dependence 
of the bare tunneling rates 7/. (hi) the transition from 
region D to G is governed by the condition: 



r (2) 
1 01 

r (2) 
1 00 



= A 2 < 



2(a + l) 
7 -2a' 



(38) 



which explains the persistence of a PDC line also for 
smaller values A (compare again Fig. [7] with Fig. [TO]) . 
The corrections introduced by a finite Tii rate do not 
change qualitatively the analysis and can be found for 
completeness in the Appendix [Bl 



C. Effect of an asymmetric coupling of the 
different orbital states to the leads 

The Franck Condon factors are still assumed to be the 
same for the source and drain tunneling and no overall 



V. EFFECT OF THE ASYMMETRIC 
COUPLING TO THE LEFT AND RIGHT LEAD 

Eventually, we consider the effect of an asymmetry in 
the coupling to the left and right lead in combination 
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to the orbital asymmetry discussed in the previous sec- 
tion. We introduce the asymmetry via the parameter 
b = Jlj/irj where I = 1, 2 represents the orbital degree 
of freedom and by convention j Sj i = for positive bias 
voltages. In Fig. [IT] we present the stability diagrams for 
a molecule coupled to vibrons with an orbital asymmetry 
a = 1/45 and two different left and right asymmetries. 
In Fig. Djja) the asymmetry parameter b = 45 while 
b = 1/45 in Fig. Qjjb). Since b is the only parameter 
of the system that breaks the left and right symmetry, 
the differential conductances in Fig. [11] can be obtained 
one from the other by a reflection of the bias. The most 
striking effect of the left and right asymmetry is, though, 
to produce the NDC lines always in the same direction 
[compare Figs.fTQlc) andQjJa)] both for positive and neg- 
ative biases. As in the previous sections, we studied ana- 
lytically the transitions lines separating the A, B, D and 
G low-bias regions of the stability diagram (see Fig. [8]). 
We could thus obtain the NDC conditions for arbitrary 
values of the asymmetry parameters a and b. The tran- 
sition line between region A and B is an NDC line for 
every electron phonon coupling A under the condition 



a > 



2 b +1 
~2b~ " 



(39) 



The NDC condition for the transition between the B and 
D region reads instead 



A < 



2(o+l) 



46' 



(40) 



while the transition between the regions D and G is gov- 
erned by the relation 



A < 



2(a+ 1)6 
l + 2(3-a)6' 



(41) 



Equations (|39|) - (|4l]) allow a partial interpretation of the 
numerical results presented in Fig. [TTJ It is in fact easy 
to demonstrate that, if b > 1, for sufficiently small values 
of a (a < 1/6) the alternating NDC and PDC pattern 
at positive biases is not modified by the left and right 
asymmetry parameter b . With the help of the same 
set of equations and the symmetry property mentioned 
above we can also analyze the negative bias transitions. 
The sequence of transitions between the regions A, B, D, 
and G in Fig. [HT b) reveals in fact a very different pat- 
tern of strong PDC transitions alternated by very weak 
NDC lines. This sequence can be obtained by the condi- 
tions expressed in Eqs. (|39|) - (|^T|) by substituting b with 
1/6 in the limit b ^> 1. Unfortunately, due to the pe- 
culiar structure of the harmonic spectrum, the number 
of states involved in the regions E and F of the stability 
diagram grows rapidly and the analytical analysis of the 
transition, even if possible, becomes very cumbersome. 
Our numerical findings are nevertheless consistent with 
the ones reported by other groups^, where the left and 
right asymmetry has been given as a necessary condition 
for achieving NDC with the same slope at both positive 
and negative biases. 
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Figure 11: (Color online). (a),(b) Stability diagrams for both 
an orbital and left and right asymmetry. All the parameters 
are the same as used in Fig. flQl c). The asymmetry with 
respect to the left and right lead is b = 45 in panel (a) while 
b = 1/45 in panel (b). 



VI. CONCLUSIONS 

In this paper we analyzed the spectrum and the 
vibration-assisted transport properties of a nanostruc- 
ture where two degenerate or quasidegenerate levels are 
coupled to several vibronic modes. Our model can cap- 
ture features of transport properties of suspended carbon 
nanotubes and of molecules with a fourfold degenerate 
electronic level coupled to many vibrational modes. The 
transport theory is based on vibron-assisted tunneling, 
mediated by vibrational modes. In order to study the dy- 
namics of the system, we apply a density matrix approach 
which starts from the Liouville equation for the total den- 
sity operator and which enables the treatment of degener- 
ate and quasidegenerate states. Despite the fact that we 
considered a fully symmetric setup, the stability diagram 
for the differential conductance shows striking negative 
differential conductance (NDC) features which hints at 
peculiar features of our nanoelectromechanical system. 
We predict that NDCs appear due to the slow channels 
in the source transition originating from spin and/or or- 
bital degeneracies and the suppression of Franck- Condon 
channels. With source-drain symmetry being preserved 
but an asymmetry between orbitally degenerate states 
being allowed, we could explain the alternating PDC and 
NDC features observed in Refs.0and[9|. Eventually, with 
the further introduction of the left and right asymmetry 
suggested in Ref. 12, we confirmed the appearance, also 
in presence of multimodes, of the NDC lines with the 
same slope for both positive and negative biases. We 
also gave an analytical interpretation of the numerical 
results in the low-bias regime. 



Acknowledgments 

We acknowledge the support of DFG under the pro- 
grams SFB 689 and GRK 1570. A. Y. acknowledges 
the support of Kohat University of Science & Technol- 
ogy, Kohat-26000, Khyber Pakhtunkhwa, Pakistan under 
the Human Resource Development Program. We thank 
Georg Begemann and Leonhard Mayrhofer for their help 



10 



in the numerical calculations. 



where m 



mm max 



min/max(m, m'). 



Appendix A: Evaluation of transition matrix 
elements of electron operator 



Appendix B: NDC and/or PDC threshold 
conditions with finite Fn 



To determine the transition rates, we calculate the ma- 
trix elements 

(r\d la \s) = e-s ^l A -l 2 J] F (A n , m n , m' n ) , (Al) 

n 

where \r) and \s) represent the eigenstates given by 
Eq. (fTD|) . The function F(A,ra,ra') determines the cou- 
pling between states with a different vibronic number of 
excitations with effective coupling A and is expressed as 

F(A, m, ra') = fe(m' - m)X m '- m + 9(m - ra') 

(-|A| 2 )^ 



:(-A*) m ~ m/ ) 



-If,". 



TTlffiax ' Z- ~^ z!(z -|- Tfiffiax ^mm)« 



(A2) 



In this appendix we give the conditions which deter- 
mines the sign of the current change in the transition 
from region B to D and D to G with finite mismatch 6a 
(Fig. |8]). We take only the first order contribution in the 
ratios Fh/Fqq and Fh/Fqi. The validity of these for- 
mulas is thus restricted to A w 1. The condition for the 
transition B to D reads 



-(2) 
01 



< 



2 + 2a, 



5a 



^(2) 
00 



14 + 9a (2) 
20(1 + a) 11 ' 



(Bl) 



while for the transition D to G one obtains 

,(2) ^ 2(o+l) (2) _4o 2 -16o-ll (2 ) 



0! 7 . 



2a 



00 



4(a+ l)(2a- 7) 



rff, (B2) 



where a = r ) a i/la2 measures in both cases the asymme- 
try between the coupling to the different orbit als. 
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